#include "sgPArpack.h"
// #include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <iostream>
#include <limits>  // std::numeric_limits
#include <memory>

#include "parpack.h"

#include "Factory.h"
#include "LinearSolver.h"
#include "SparseBlas.h"
#include "debug_c.h"  // debug arpack.
#include "stat_c.h"   // arpack statistics.
#include <complex.h>  // creal, cimag.

using namespace SG::Algebra;

// Alias to avoid the tiresome long names.
typedef std::vector<double> DoubleArray;
typedef std::vector<a_int>  IntArray;

#define FREE_MEMORY(A) \
    do                 \
    {                  \
        if (A)         \
        {              \
            free (A);  \
        }              \
    } while (0)

namespace SG::Algebra
{
    template <typename TDataType>
    struct PArpackPrivate
    {
        MPI_Fint                            comm;  // MPI  Communicator for the processor grid
        int                                 info;
        int                                 n;  // Dimension of the eigenproblem.
        typename PArpack<TDataType>::Option opt;
    };

    template <typename TDataType>
    std::shared_ptr<TLinearSolver<TDataType>> getLinearSolver (const std::string& name)
    {
        std::shared_ptr<TLinearSolver<TDataType>> ls;

        auto factory = SG::Base::Factory::GetInstance ();
        ls.reset (( TLinearSolver<TDataType>* ) (factory->Produce (name)));

        return ls;
    }

    char const* convert_to_char (const arpack::BMat& option)
    {
        return option == arpack::BMat::Identity ? "I" : "G";
    }

    char const* convert_to_char (const arpack::Which& option)
    {
        using namespace arpack;

        switch (option)
        {
        case Which::LargestAlgebraic: {
            return "LA";
            break;
        }
        case Which::SmallestAlgebraic: {
            return "SA";
            break;
        }
        case Which::LargestMagnitude: {
            return "LM";
            break;
        }
        case Which::SmallestMagnitude: {
            return "SM";
            break;
        }
        case Which::LargestReal: {
            return "LR";
            break;
        }
        case Which::SmallestReal: {
            return "SR";
            break;
        }
        case Which::LargestImaginary: {
            return "LI";
            break;
        }
        case Which::SmallestImaginary: {
            return "SI";
            break;
        }
            // case Which::both_ends: {
            //     return "BE";
            //     break;
            // }
        }
        return "LM";
    }

    /**
    c  HOWMNY  Character*1  (INPUT)
    c          Specifies the form of the basis for the invariant subspace
    c          corresponding to the converged Ritz values that is to be computed.
    c
    c          = 'A': Compute NEV Ritz vectors;
    c          = 'P': Compute NEV Schur vectors;
    c          = 'S': compute some of the Ritz vectors, specified
    c                 by the logical array SELECT.
     */
    enum class HowMany : int
    {
        /// 'A' Compute NEV Ritz vectors
        ritz_vectors,
        /// 'P' Compute NEV Schur vectors;
        schur_vectors,
        /// 'S' compute some of the Ritz vectors, specified by the logical array
        /// SELECT.
        ritz_specified
    };

    char const* convert_to_char (const HowMany& option)
    {
        switch (option)
        {
        case HowMany::ritz_vectors: {
            return "A";
            break;
        }
        case HowMany::schur_vectors: {
            return "P";
            break;
        }
        case HowMany::ritz_specified: {
            return "S";
            break;
        }
        }
        return "A";
    }

}  // namespace SG::Algebra

//=============================================================================
// sgPArpack

template <typename TDataType>
PArpack<TDataType>::PArpack ()
{
    imp       = new PArpackPrivate<TDataType> ();
    imp->comm = MPI_Comm_c2f (MPI_COMM_WORLD);
}

template <typename TDataType>
PArpack<TDataType>::~PArpack ()
{
    if (imp)
    {
        delete imp;
    }
}

template <typename TDataType>
typename PArpack<TDataType>::Option PArpack<TDataType>::getDefaultOption ()
{
    // GoF's Prototype Pattern
    // create new objects by copying this prototype.
    static PArpack<TDataType>::Option opt;
    // opt.which = "LA";
    // opt.nev   = 1;
    // opt.sigma = 0.0;

    return opt;
}

template <typename TDataType>
void PArpack<TDataType>::setOption (const Option& opt)
{
    imp->opt = opt;
}

template <typename TDataType>
Info_t PArpack<TDataType>::info () const
{
    return ( Info_t ) imp->info;
}

//=============================================================================
// Explicit instantiations to hide the underlying implementation.
namespace SG
{
    namespace Algebra
    {
        template class ALGEBRA_EXPORT PArpack<Real_t>;
        template class ALGEBRA_EXPORT PArpack<Complex_t>;
    }  // namespace Algebra
}  // namespace SG

//=============================================================================
// Explicit specilizations

//------------------------------------------------------------------
// Real_t

namespace  // Anonymous
{
    /*
     * Re-formed from arpack-ng/PARPACK/TESTS/MPI/icb_parpack_c.c
     *                dealii/include/deal.II/lac/parpack_solver.h
     *
     * NOTE: If the linear operator "OP" is real and symmetric
     * with respect to the real positive semi-definite symmetric matrix B,
     * i.e. B*OP = (OP`)*B, then subroutine dsaupd  should be used instead.
     *
     * Several key notes are shown as the followings,
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pdsaupd.f
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pdseupd.f
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pdnaupd.f
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pdneupd.f
     *
     */

    /*
     * Mode 2:  A*x = lambda*M*x, A symmetric, M symmetric positive definite
     *      ===> OP = inv[M]*A  and  B = M.
     *       ===> (If M can be factored see remark 3 below)
     */
    int eig2_d (const PArpackPrivate<double>* imp,
                const SpMat_t&                K,
                const SpMat_t&                M,
                std::vector<Real_t>&          eigenvalues,
                std::vector<Vec_t>&           eigenvectors)
    {
        int ret = 0;

        double* resid  = nullptr;
        double* V      = nullptr;
        double* z      = nullptr;
        double* d      = nullptr;
        double* workd  = nullptr;
        double* workl  = nullptr;
        a_int*  select = nullptr;

        MPI_Fint    comm = imp->comm;  // MPI  Communicator for the processor grid

        const a_int N = K.cols ();
        const a_int N_local   = K.local_rows ();
        const a_int row_start = K.row_start ();
        const a_int row_end   = K.row_end ();

        const a_int nev = imp->opt.nev;
        // how many Lanczos vectors are generated at each iteration
        // see octave/liboctave/numeric/eigs-base.cc
        // see deal.ii/include/deal.II/lac/arpack_solver.h
        const a_int ncv    = std::min (std::max (2 * nev + 1, 20), N);
        const a_int ldv    = N_local;
        const a_int ldz    = N_local;
        const a_int lworkl = ncv * (ncv + 8);
        const a_int rvec   = 1;         // need eigenvectors

        const double tol   = 0.000001;  // small tol => more stable checks after EV computation.
        const double sigma = 0;         // not referenced in this mode

        // Maximum number of iterations.
        // see octave/libinterp/__eigs__.cc
        const a_int maxit = 300;

        resid = ( double* ) malloc (N_local * sizeof (double));
        V     = ( double* ) malloc (ldv * ncv * sizeof (double));
        z     = ( double* ) malloc (ldz * nev * sizeof (double));
        d     = ( double* ) malloc (nev * sizeof (double));
        workd = ( double* ) malloc (3 * N_local * sizeof (double));
        workl = ( double* ) malloc (lworkl * sizeof (double));

        select = ( a_int* ) malloc (ncv * sizeof (a_int));  // since HOWMNY = 'A', only used as workspace here

        a_int iparam[11], ipntr[11];
        iparam[0] = 1;                                      // ishift
        // iparam[2] = 10 * N;                                 // on input: maxit; on output: actual iteration
        iparam[2] = maxit;
        iparam[3] = 1;  // NB, only 1 allowed
        iparam[6] = 2;  // mode

        // char bmat[] = "G";                                  // generalized eigenvalue problem A*x = lambda*B*x
        // const char* bmat = "G";  // generalized eigenvalue problem A*x = lambda*B*x
        const char* bmat = convert_to_char (imp->opt.bmat);  // generalized eigenvalue problem A*x = lambda*B*x
        /*
c\Remarks
c  1. The converged Ritz values are always returned in ascending
c     algebraic order.  The computed Ritz values are approximate
c     eigenvalues of OP.  The selection of WHICH should be made
c     with this in mind when Mode = 3,4,5.  After convergence,
c     approximate eigenvalues of the original problem may be obtained
c     with the ARPACK subroutine dseupd .
        */
        // char which[]  = "LM";                               // compute the NEV largest (in magnitude) eigenvalues
        const char* which = imp->opt.which.c_str ();
        // char        howmny[] = "A";
        const char* howmny = "A";

        // consider sigma = 0 for now.
#if defined(SG_USE_MKL)
        auto ls = getLinearSolver<double> ("CPardiso<Real_t>");
#else
        auto ls = getLinearSolver<double> ("TPetscKsp<Real_t>");
#endif
        ls->compute (M);

        a_int info = 0, ido = 0;
        do
        {
            pdsaupd_c (comm, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);

            // solve Y = OP * X
            // dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
            if (ido == -1)
            {
                /*
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t kx = K * x;

                Vec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 1)
            {
                /*
c          IDO =  1: compute  Y = OP * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    In mode 3,4 and 5, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t kx = K * x;

                Vec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 2)
            {
                /*
c          IDO =  2: compute  Y = B * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t mx = M * x;

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = mx[i];
                }
            }
            else if (ido != 99)
            {
                printf ("Error in dsaupd: unexpected ido %d\n", ido);
                ret = 1;
                goto EXIT;
            }

        } while (ido != 99);

        printf ("dsaupd: number of Arnoldi update iterations %d/%d\n", iparam[2], maxit);
        printf ("dsaupd: NUMOP %d, NUMOPB %d, NUMREO %d\n", iparam[8], iparam[9], iparam[10]);

        // check info and number of ev found by arpack.
        if (info < 0 || iparam[4] < nev)
        {
            printf ("Error in dsaupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
            ret = 1;
            goto EXIT;
        }

        pdseupd_c (comm,
                   rvec,
                   howmny,
                   select,
                   d,
                   z,
                   ldz,
                   sigma,
                   bmat,
                   N_local,
                   which,
                   nev,
                   tol,
                   resid,
                   ncv,
                   V,
                   ldv,
                   iparam,
                   ipntr,
                   workd,
                   workl,
                   lworkl,
                   &info);
        if (info < 0)
        {
            printf ("Error in dsaupd: info %d\n", info);
            ret = 1;
            goto EXIT;
        }

        // store the eigenvectors and eigenvalues
        eigenvalues.resize (nev);
        eigenvectors.resize (nev);

        for (int j = 0; j < nev; ++j)
        {
            eigenvectors[j].resize (N);

            for (int i = row_start; i < row_end; ++i)
            {
                eigenvectors[j].set(i, z[j * N_local + i - row_start]);
            }

            eigenvectors[j].assemble ();

            eigenvalues[j] = d[j];
        }

    EXIT:
        FREE_MEMORY (resid);
        FREE_MEMORY (V);
        FREE_MEMORY (z);
        FREE_MEMORY (d);
        FREE_MEMORY (workd);
        FREE_MEMORY (workl);
        FREE_MEMORY (select);

        return ret;
    }

    /*
     * Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite.
     *          ===> OP = (inv[K - sigma*M])*M  and  B = M.
     *          ===> Shift-and-Invert mode
     */
    int eig3_d (const PArpackPrivate<double>* imp,
                const SpMat_t&                K,
                const SpMat_t&                M,
                std::vector<Real_t>&          eigenvalues,
                std::vector<Vec_t>&           eigenvectors)
    {
        int ret = 0;

        double* resid  = nullptr;
        double* V      = nullptr;
        double* z      = nullptr;
        double* d      = nullptr;
        double* workd  = nullptr;
        double* workl  = nullptr;
        a_int*  select = nullptr;

        MPI_Fint    comm = imp->comm;  // MPI  Communicator for the processor grid

        const a_int N = K.cols ();
        const a_int N_local   = K.local_rows ();
        const a_int row_start = K.row_start ();
        const a_int row_end   = K.row_end ();

        const a_int nev = imp->opt.nev;
        // how many Lanczos vectors are generated at each iteration
        // see octave/liboctave/numeric/eigs-base.cc
        // see deal.ii/include/deal.II/lac/arpack_solver.h
        const a_int ncv    = std::min (std::max (2 * nev + 1, 20), N);
        const a_int ldv    = N_local;
        const a_int ldz    = N_local;
        const a_int lworkl = ncv * (ncv + 8);
        const a_int rvec   = 1;       // need eigenvectors

        const double tol = 0.000001;  // small tol => more stable checks after EV computation.
        // const double sigma = 0;         // represents the shift
        const double sigma = imp->opt.sigma;  // represents the shift

        // Maximum number of iterations.
        // see octave/libinterp/__eigs__.cc
        const a_int maxit = 300;

        resid = ( double* ) malloc (N_local * sizeof (double));
        V     = ( double* ) malloc (ldv * ncv * sizeof (double));
        z     = ( double* ) malloc (ldz * nev * sizeof (double));
        d     = ( double* ) malloc (nev * sizeof (double));
        workd = ( double* ) malloc (3 * N_local * sizeof (double));
        workl = ( double* ) malloc (lworkl * sizeof (double));

        select = ( a_int* ) malloc (ncv * sizeof (a_int));  // since HOWMNY = 'A', only used as workspace here

        a_int iparam[11], ipntr[11];
        iparam[0] = 1;                                      // ishift
        // iparam[2] = 10 * N;                                 // on input: maxit; on output: actual iteration
        iparam[2] = maxit;
        iparam[3] = 1;  // NB, only 1 allowed
        iparam[6] = 3;  // mode

        // char bmat[] = "G";                                  // generalized eigenvalue problem A*x = lambda*B*x
        // const char* bmat = "G";  // generalized eigenvalue problem A*x = lambda*B*x
        const char* bmat = convert_to_char (imp->opt.bmat);  // generalized eigenvalue problem A*x = lambda*B*x
        /*
c\Remarks
c  1. The converged Ritz values are always returned in ascending
c     algebraic order.  The computed Ritz values are approximate
c     eigenvalues of OP.  The selection of WHICH should be made
c     with this in mind when Mode = 3,4,5.  After convergence,
c     approximate eigenvalues of the original problem may be obtained
c     with the ARPACK subroutine dseupd .
        */
        // char which[]  = "LM";                               // compute the NEV largest (in magnitude) eigenvalues
        const char* which = imp->opt.which.c_str ();
        // char        howmny[] = "A";
        const char* howmny = "A";

#if defined(SG_USE_MKL)
        auto ls = getLinearSolver<double> ("CPardiso<Real_t>");
#else
        auto ls = getLinearSolver<double> ("TPetscKsp<Real_t>");
#endif

        if (sigma == 0.0)  // consider sigma = 0
        {
            ls->compute (K);
        }
        else  // consider sigma != 0
        {
            SpMat_t Ksigma = K - sigma * M;
            ls->compute (Ksigma);
        }

        a_int info = 0, ido = 0;
        do
        {
            pdsaupd_c (comm, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);

            // solve Y = OP * X
            // dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
            if (ido == -1)
            {
                /*
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t mx = M * x;

                Vec_t y;
                ls->solve (mx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 1)
            {
                /*
c          IDO =  1: compute  Y = OP * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    In mode 3,4 and 5, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.
                */
                Vec_t mx (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    mx.set (i, workd[ipntr[2] - 1 + i - row_start]);
                }
                mx.assemble ();

                Vec_t y;
                ls->solve (mx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 2)
            {
                /*
c          IDO =  2: compute  Y = B * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t mx = M * x;

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = mx[i];
                }
            }
            else if (ido != 99)
            {
                printf ("Error in dsaupd: unexpected ido %d\n", ido);
                ret = 1;
                goto EXIT;
            }

        } while (ido != 99);

        printf ("dsaupd: number of Arnoldi update iterations %d/%d\n", iparam[2], maxit);
        printf ("dsaupd: NUMOP %d, NUMOPB %d, NUMREO %d\n", iparam[8], iparam[9], iparam[10]);

        // check info and number of ev found by arpack.
        if (info < 0 || iparam[4] < nev)
        {
            printf ("Error in dsaupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
            ret = 1;
            goto EXIT;
        }

        pdseupd_c (comm,
                   rvec,
                   howmny,
                   select,
                   d,
                   z,
                   ldz,
                   sigma,
                   bmat,
                   N_local,
                   which,
                   nev,
                   tol,
                   resid,
                   ncv,
                   V,
                   ldv,
                   iparam,
                   ipntr,
                   workd,
                   workl,
                   lworkl,
                   &info);
        if (info < 0)
        {
            printf ("Error in dsaupd: info %d\n", info);
            ret = 1;
            goto EXIT;
        }

        // store the eigenvectors and eigenvalues
        eigenvalues.resize (nev);
        eigenvectors.resize (nev);

        for (int j = 0; j < nev; ++j)
        {
            eigenvectors[j].resize (N);

            for (int i = row_start; i < row_end; ++i)
            {
                eigenvectors[j].set (i, z[j * N_local + i - row_start]);
            }

            eigenvectors[j].assemble ();

            eigenvalues[j] = d[j];
        }

    EXIT:
        FREE_MEMORY (resid);
        FREE_MEMORY (V);
        FREE_MEMORY (z);
        FREE_MEMORY (d);
        FREE_MEMORY (workd);
        FREE_MEMORY (workl);
        FREE_MEMORY (select);

        return ret;
    }

    /*
     * Mimic the Buckling mode using Mode 3.
     *
     * ARPACK Mode 4 seems un-stable to check out the tageted eigenvalues. What's more,
     * the shift value 'sigma' must be given non-zero value, namely, sigma != 0.
     *
     * @note the resultant eigenvalues will be presented in \f$ 1/lamda \f$ ascending algebraic order.
     *
     * Mode 4:  K*x = lambda*KG*x, K symmetric positive semi-definite,
     *          KG symmetric indefinite
     *          ===> OP = (inv[KG - sigma*K])*K  and  B = K.
     */
    int eig4_d (const PArpackPrivate<double>* imp,
                const SpMat_t&                K,
                const SpMat_t&                KG,
                std::vector<Real_t>&          eigenvalues,
                std::vector<Vec_t>&           eigenvectors)
    {
        int ret = 0;

        double* resid  = nullptr;
        double* V      = nullptr;
        double* z      = nullptr;
        double* d      = nullptr;
        double* workd  = nullptr;
        double* workl  = nullptr;
        a_int*  select = nullptr;

        MPI_Fint comm = imp->comm;  // MPI  Communicator for the processor grid

        const a_int N = K.cols ();
        const a_int N_local   = K.local_rows ();
        const a_int row_start = K.row_start ();
        const a_int row_end   = K.row_end ();

        const a_int nev = imp->opt.nev;
        // how many Lanczos vectors are generated at each iteration
        // see octave/liboctave/numeric/eigs-base.cc
        // see deal.ii/include/deal.II/lac/arpack_solver.h
        const a_int ncv    = std::min (std::max (2 * nev + 1, 20), N);
        const a_int ldv    = N_local;
        const a_int ldz    = N_local;
        const a_int lworkl = ncv * (ncv + 8);
        const a_int rvec   = 1;       // need eigenvectors

        const double tol = 0.000001;  // small tol => more stable checks after EV computation.
        // const double sigma = 0;         // represents the shift
        const double sigma = imp->opt.sigma;  // represents the shift

        // Maximum number of iterations.
        // see octave/libinterp/__eigs__.cc
        const a_int maxit = 300;

        resid = ( double* ) malloc (N_local * sizeof (double));
        V     = ( double* ) malloc (ldv * ncv * sizeof (double));
        z     = ( double* ) malloc (ldz * nev * sizeof (double));
        d     = ( double* ) malloc (nev * sizeof (double));
        workd = ( double* ) malloc (3 * N_local * sizeof (double));
        workl = ( double* ) malloc (lworkl * sizeof (double));

        select = ( a_int* ) malloc (ncv * sizeof (a_int));  // since HOWMNY = 'A', only used as workspace here

        a_int iparam[11], ipntr[11];
        iparam[0] = 1;                                      // ishift
        // iparam[2] = 10 * N;                                 // on input: maxit; on output: actual iteration
        iparam[2] = maxit;
        iparam[3] = 1;  // NB, only 1 allowed
        // iparam[6] = 4;  // mode
        iparam[6] = 3;  // mode: use mode 3 to mimic buckling mode

        // char bmat[] = "G";                                  // generalized eigenvalue problem A*x = lambda*B*x
        // const char* bmat = "G";  // generalized eigenvalue problem A*x = lambda*B*x
        const char* bmat = convert_to_char (imp->opt.bmat);  // generalized eigenvalue problem A*x = lambda*B*x
        /*
c\Remarks
c  1. The converged Ritz values are always returned in ascending
c     algebraic order.  The computed Ritz values are approximate
c     eigenvalues of OP.  The selection of WHICH should be made
c     with this in mind when Mode = 3,4,5.  After convergence,
c     approximate eigenvalues of the original problem may be obtained
c     with the ARPACK subroutine dseupd .
        */
        // char which[]  = "LM";                               // compute the NEV largest (in magnitude) eigenvalues
        const char* which = imp->opt.which.c_str ();
        // char        howmny[] = "A";
        const char* howmny = "A";

#if defined(SG_USE_MKL)
        auto ls = getLinearSolver<double> ("CPardiso<Real_t>");
#else
        auto ls = getLinearSolver<double> ("TPetscKsp<Real_t>");
#endif

        // if (sigma == 0.0)  // consider sigma = 0
        // {
        //     ls->compute (K);
        // }
        // else  // consider sigma != 0
        // {
        //     SpMat_t Ksigma = K - sigma * KG;
        //     ls->compute (Ksigma);
        // }

        // SpMat_t Ksigma = K - sigma * KG;
        SpMat_t Ksigma = KG - sigma * K;
        ls->compute (Ksigma);

        a_int info = 0, ido = 0;
        do
        {
            pdsaupd_c (comm, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);

            // solve Y = OP * X
            // dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
            if (ido == -1)
            {
                /*
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t kx = K * x;

                Vec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 1)
            {
                /*
c          IDO =  1: compute  Y = OP * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    In mode 3,4 and 5, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.
                */
                Vec_t kx (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    kx.set(i, workd[ipntr[2] - 1 + i - row_start]);
                }
                kx.assemble ();

                Vec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = y[i];
                }
            }
            else if (ido == 2)
            {
                /*
c          IDO =  2: compute  Y = B * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
                */
                Vec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                Vec_t kx = K * x;

                for (int i = row_start; i < row_end; ++i)
                {
                    workd[ipntr[1] - 1 + i - row_start] = kx[i];
                }
            }
            else if (ido != 99)
            {
                printf ("Error in dsaupd: unexpected ido %d\n", ido);
                ret = 1;
                goto EXIT;
            }

        } while (ido != 99);

        printf ("dsaupd: number of Arnoldi update iterations %d/%d\n", iparam[2], maxit);
        printf ("dsaupd: NUMOP %d, NUMOPB %d, NUMREO %d\n", iparam[8], iparam[9], iparam[10]);

        // check info and number of ev found by arpack.
        if (info < 0 || iparam[4] < nev)
        {
            printf ("Error in dsaupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
            ret = 1;
            goto EXIT;
        }

        pdseupd_c (comm,
                   rvec,
                   howmny,
                   select,
                   d,
                   z,
                   ldz,
                   sigma,
                   bmat,
                   N_local,
                   which,
                   nev,
                   tol,
                   resid,
                   ncv,
                   V,
                   ldv,
                   iparam,
                   ipntr,
                   workd,
                   workl,
                   lworkl,
                   &info);
        if (info < 0)
        {
            printf ("Error in dsaupd: info %d\n", info);
            ret = 1;
            goto EXIT;
        }

        // store the eigenvectors and eigenvalues
        eigenvalues.resize (nev);
        eigenvectors.resize (nev);

        for (int j = 0; j < nev; ++j)
        {
            eigenvectors[j].resize (N);

            for (int i = row_start; i < row_end; ++i)
            {
                eigenvectors[j].set (i, z[j * N_local + i - row_start]);
            }

            eigenvectors[j].assemble ();

            eigenvalues[j] = d[j];
        }

        // std::cout << "E = " << std::endl;
        // std::cout << E << std::endl;

    EXIT:
        FREE_MEMORY (resid);
        FREE_MEMORY (V);
        FREE_MEMORY (z);
        FREE_MEMORY (d);
        FREE_MEMORY (workd);
        FREE_MEMORY (workl);
        FREE_MEMORY (select);

        return ret;
    }

}  // namespace

template <>
void ALGEBRA_EXPORT PArpack<Real_t>::solve (const SparseMatrixType&  A,
                                            const SparseMatrixType&  B,
                                            std::vector<Real_t>&     eigenvalues,
                                            std::vector<VectorType>& eigenvectors)
{
    imp->comm = MPI_Comm_c2f (MPI_COMM_WORLD);
    imp->info = 0;
    imp->n    = A.rows ();
    // imp->nev   = 1;
    // imp->opt.sigma = 0.0;

    // checking the options...
    // see octave/liboctave/numeric/eigs-base.cc
    if (imp->opt.nev < 1 || imp->opt.nev > imp->n - 1)
    {
        printf ("Invalid number of eigenvalues to extract (must be 0 < k < n) \n");
        return;
    }

    // [eigenvectors, eigenvalues]
    // eigenvectors: E(:. 0:nev)
    // eigenvalues: E(0:nev, nev)
    // Mat_t E (imp->n, imp->opt.nev + 1);
    if (imp->opt.mode == arpack::Mode::MODE_2)
    {
        imp->info = eig2_d (imp, A, B, eigenvalues, eigenvectors);
    }
    else if (imp->opt.mode == arpack::Mode::MODE_3)
    {
        imp->info = eig3_d (imp, A, B, eigenvalues, eigenvectors);
    }
    else if (imp->opt.mode == arpack::Mode::MODE_4)
    {
        imp->info = eig4_d (imp, A, B, eigenvalues, eigenvectors);
    }
}

//------------------------------------------------------------------
// Complex_t

//
// some helpful notes are listed as followings,
// 1. complex processing in GCC
// ISO C99 supports complex floating data types, and as an extension GCC supports them in C90 mode and in C++. GCC also supports complex
// integer data types which are not part of ISO C99. You can declare complex types using the keyword _Complex. As an extension, the
// older GNU keyword __complex__ is also supported. For example, ‘_Complex double x;’ declares x as a variable whose real part and
// imaginary part are both of type double. ‘_Complex short int y;’ declares y to have real and imaginary parts of type short int; this
// is not likely to be useful, but it shows that the set of complex types is complete.

// To write a constant with a complex data type, use the suffix ‘i’ or ‘j’ (either one; they are equivalent).For example,
// 2.5fi has type _Complex float and 3i has type _Complex int.Such a constant always has a pure imaginary value,
// but you can form any complex value you like by adding one to a real constant.This is a GNU extension;
// if you have an ISO C99 conforming C library (such as the GNU C Library), and want to construct complex constants of floating type,
// you should include<complex.h>and use the macros I or _Complex_I instead.
//
// see https://gcc.gnu.org/onlinedocs/gcc-7.5.0/gcc/Complex.html
//
// 2. the 'const' type qualifier for the character array required by the 'znaupd/zneupd' routines
//
// On Linux, the character arrays without of that 'const' type qualifier, seem to be modified after the 'znaupd/zneupd' routines return.
// Make sure that 'bamt/which/howmny' strings are qualified with the 'const'.
//
// the debug clues with VSCode on Linux: eig3_z=>znaupd_c=>icbazn.F90=>znaupd.f=>icbazn.F90=>znaupd_c=>eig3_z
//
//

namespace
{
    /*
     * Re-formed from arpack-ng/TESTS/icb_arpck_c.cc
     *
     * Several key notes are shown as the followings,
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pznaupd.f
     * see https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pzneupd.f
     *
     */

    /**
     * Mode 2:  A*x = lambda*M*x, M hermitian positive definite
     *            ===> OP = inv[M]*A  and  B = M.
     *            ===> (If M can be factored see remark 3 below)
     */
    int eig2_z (const PArpackPrivate<Complex_t>* imp,
                const ComplexSpMat_t&            K,
                const ComplexSpMat_t&            M,
                std::vector<Complex_t>&          eigenvalues,
                std::vector<ComplexVec_t>&       eigenvectors)
    {
        int ret = 0;

        a_dcomplex* resid  = nullptr;
        a_dcomplex* V      = nullptr;
        a_dcomplex* z      = nullptr;
        a_dcomplex* d      = nullptr;
        a_dcomplex* workd  = nullptr;
        a_dcomplex* workl  = nullptr;
        a_dcomplex* workev = nullptr;
        double*     rwork  = nullptr;
        a_int*      select = nullptr;

        MPI_Fint comm = imp->comm;  // MPI  Communicator for the processor grid

        const a_int N         = K.cols ();
        const a_int N_local   = K.local_rows ();
        const a_int row_start = K.row_start ();
        const a_int row_end   = K.row_end ();

        const a_int nev = imp->opt.nev;
        // how many Lanczos vectors are generated at each iteration
        // see octave/liboctave/numeric/eigs-base.cc
        // see deal.ii/include/deal.II/lac/arpack_solver.h
        const a_int ncv    = std::min (std::max (2 * nev + 1, 20), N);
        const a_int ldv    = N_local;
        const a_int ldz    = N_local;
        const a_int lworkl = ncv * (3 * ncv + 5);
        const a_int rvec   = 1;       // need eigenvectors

        const double tol = 0.000001;  // small tol => more stable checks after EV computation.
        // const a_dcomplex sigma = CMPLX (0., 0.);  // not referenced in this mode
        const a_dcomplex sigma = 0.0;
        memcpy (( void* ) &sigma, ( void* ) &imp->opt.sigma, sizeof (a_dcomplex));

        // Maximum number of iterations.
        // see octave/libinterp/__eigs__.cc
        const a_int maxit = 300;

        resid  = ( a_dcomplex* ) malloc (N_local * sizeof (a_dcomplex));
        V      = ( a_dcomplex* ) malloc (ldv * ncv * sizeof (a_dcomplex));
        z      = ( a_dcomplex* ) malloc (ldz * nev * sizeof (a_dcomplex));
        d      = ( a_dcomplex* ) malloc (nev * sizeof (a_dcomplex));
        workd  = ( a_dcomplex* ) malloc (3 * N_local * sizeof (a_dcomplex));
        workl  = ( a_dcomplex* ) malloc (lworkl * sizeof (a_dcomplex));
        workev = ( a_dcomplex* ) malloc (2 * ncv * sizeof (a_dcomplex));

        rwork  = ( double* ) malloc (ncv * sizeof (double));
        select = ( a_int* ) malloc (ncv * sizeof (a_int));  // since HOWMNY = 'A', only used as workspace here

        a_int iparam[11], ipntr[11];
        iparam[0] = 1;                                      // ishift
        // iparam[2] = 10 * N;                                 // on input: maxit; on output: actual iteration
        iparam[2] = maxit;
        iparam[3] = 1;  // NB, only 1 allowed
        iparam[6] = 2;  // mode

        // char bmat[] = "G";                                  // generalized eigenvalue problem A*x = lambda*B*x
        // const char* bmat = "G";  // generalized eigenvalue problem A*x = lambda*B*x
        const char* bmat = convert_to_char (imp->opt.bmat);  // generalized eigenvalue problem A*x = lambda*B*x

        /*
c\Remarks
c  1. The converged Ritz values are always returned in ascending
c     algebraic order.  The computed Ritz values are approximate
c     eigenvalues of OP.  The selection of WHICH should be made
c     with this in mind when Mode = 3,4,5.  After convergence,
c     approximate eigenvalues of the original problem may be obtained
c     with the ARPACK subroutine dseupd .
        */
        // char which[]  = "LM";                               // compute the NEV largest (in magnitude) eigenvalues
        // const char* which    = imp->opt.which.c_str ();
        // char        howmny[] = "A";
        const char* which  = imp->opt.which.c_str ();
        const char* howmny = "A";

        // consider sigma = 0 for now.
#if defined(SG_USE_MKL)
        auto ls = getLinearSolver<Complex_t> ("CPardiso<Complex_t>");
#else
        auto ls = getLinearSolver<Complex_t> ("TPetscKsp<Complex_t>");
#endif

        ls->compute (M);

        a_int info = 0, ido = 0;

        // {
        //     // If INFO .EQ. 0, a randomly initial residual vector is used.
        //     // If INFO .NE. 0, RESID contains the initial residual vector, possibly from a previous run.
        //     info     = 1;
        //     auto eps = std::numeric_limits<double>::epsilon ();
        //     for (a_int i = 0; i < N; ++ i)
        //         resid[i] = eps;
        // }

        do
        {
            pznaupd_c (comm, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork, &info);
            // znaupd_c (&ido, convert_to_char(bmat), N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork,
            // &info);

            // solve Y = OP * X
            // dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
            if (ido == -1)
            {
                /*
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.
                */
                ComplexVec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, ( Complex_t ) workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                ComplexVec_t kx = K * x;

                ComplexVec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v                         = y[i];
                    memcpy (workd + ipntr[1] - 1 + i -row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido == 1)
            {
                /*
c          IDO =  1: compute  Y = OP * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    In mode 3,4 and 5, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.
                */
                ComplexVec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set(i, ( Complex_t ) workd[ipntr[0] - 1 + i -row_start]);
                }
                x.assemble ();

                ComplexVec_t kx = K * x;

                ComplexVec_t y;
                ls->solve (kx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v = y[i];
                    memcpy (workd + ipntr[1] - 1 + i - row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido == 2)
            {
                /*
c          IDO =  2: compute  Y = B * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
                */
                ComplexVec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, ( Complex_t ) workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                ComplexVec_t kx = K * x;

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v                         = kx[i];
                    memcpy (workd + ipntr[1] - 1 + i - row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido != 99)
            {
                printf ("Error in znaupd: unexpected ido %d\n", ido);
                ret = 1;
                goto EXIT;
            }

        } while (ido != 99);

        printf ("znaupd: number of Arnoldi update iterations %d/%d\n", iparam[2], maxit);
        printf ("znaupd: NUMOP %d, NUMOPB %d, NUMREO %d\n", iparam[8], iparam[9], iparam[10]);

        // check info and number of ev found by arpack.
        if (info < 0 || iparam[4] < nev)
        {
            printf ("Error in znaupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
            ret = 1;
            goto EXIT;
        }

        pzneupd_c (comm,
                   rvec,
                   howmny,
                   select,
                   d,
                   z,
                   ldz,
                   sigma,
                   workev,
                   bmat,
                   N_local,
                   which,
                   nev,
                   tol,
                   resid,
                   ncv,
                   V,
                   ldv,
                   iparam,
                   ipntr,
                   workd,
                   workl,
                   lworkl,
                   rwork,
                   &info);

        if (info < 0)
        {
            printf ("Error in znaupd: info %d\n", info);
            ret = 1;
            goto EXIT;
        }

        // store the eigenvectors and eigenvalues
        eigenvalues.resize (nev);
        eigenvectors.resize (nev);

        for (int j = 0; j < nev; ++j)
        {
            eigenvectors[j].resize (N);

            for (int i = row_start; i < row_end; ++i)
            {
                eigenvectors[j].set(i, (Complex_t) z[j * N_local + i - row_start]);
            }

            eigenvectors[j].assemble ();

            eigenvalues[j] = ( Complex_t ) d[j];
        }

    EXIT:
        FREE_MEMORY (resid);
        FREE_MEMORY (V);
        FREE_MEMORY (z);
        FREE_MEMORY (d);
        FREE_MEMORY (workd);
        FREE_MEMORY (workl);
        FREE_MEMORY (select);

        return ret;
    }

    /*
     *   Mode 3:  A*x = lambda*M*x, M hermitian semi-definite
     *            ===> OP =  inv[A - sigma*M]*M   and  B = M.
     *            ===> shift-and-invert mode
     *            If OP*x = amu*x, then lambda = sigma + 1/amu.
     */
    int eig3_z (const PArpackPrivate<Complex_t>* imp,
                const ComplexSpMat_t&            K,
                const ComplexSpMat_t&            M,
                std::vector<Complex_t>&          eigenvalues,
                std::vector<ComplexVec_t>&       eigenvectors)
    {
        int ret = 0;

        a_dcomplex* resid  = nullptr;
        a_dcomplex* V      = nullptr;
        a_dcomplex* z      = nullptr;
        a_dcomplex* d      = nullptr;
        a_dcomplex* workd  = nullptr;
        a_dcomplex* workl  = nullptr;
        a_dcomplex* workev = nullptr;
        double*     rwork  = nullptr;
        a_int*      select = nullptr;

        MPI_Fint comm = imp->comm;  // MPI  Communicator for the processor grid

        const a_int N         = K.cols ();
        const a_int N_local   = K.local_rows ();
        const a_int row_start = K.row_start ();
        const a_int row_end   = K.row_end ();

        const a_int nev = imp->opt.nev;
        // how many Lanczos vectors are generated at each iteration
        // see octave/liboctave/numeric/eigs-base.cc
        // see deal.ii/include/deal.II/lac/arpack_solver.h
        const a_int ncv    = std::min (std::max (2 * nev + 1, 20), N);
        const a_int ldv    = N_local;
        const a_int ldz    = N_local;
        const a_int lworkl = ncv * (3 * ncv + 5);
        const a_int rvec   = 1;       // need eigenvectors

        const double tol = 0.000001;  // small tol => more stable checks after EV computation.
        /*
         * Manipulate the complex as the Complex number arithmetic,
         * @code
         * const a_dcomplex sigma = CMPLX (imp->opt.sigma.real(), imp->opt.sigma.imag());  // represents the shift
         * @endcode
         *
         * For Linux build, the gcc of the remote pipeline complains the following error,
         * error: ‘__builtin_complex’ was not declared in this scope; did you mean ‘__builtin_toupper’?
         */
        const a_dcomplex sigma = 0.0;  // represents the shift
        memcpy (( void* ) &sigma, ( void* ) &imp->opt.sigma, sizeof (a_dcomplex));

        // Maximum number of iterations.
        // see octave/libinterp/__eigs__.cc
        const a_int maxit = 300;

        resid  = ( a_dcomplex* ) malloc (N_local * sizeof (a_dcomplex));
        V      = ( a_dcomplex* ) malloc (ldv * ncv * sizeof (a_dcomplex));
        z      = ( a_dcomplex* ) malloc (ldz * nev * sizeof (a_dcomplex));
        d      = ( a_dcomplex* ) malloc (nev * sizeof (a_dcomplex));
        workd  = ( a_dcomplex* ) malloc (3 * N_local * sizeof (a_dcomplex));
        workl  = ( a_dcomplex* ) malloc (lworkl * sizeof (a_dcomplex));
        workev = ( a_dcomplex* ) malloc (2 * ncv * sizeof (a_dcomplex));

        rwork  = ( double* ) malloc (ncv * sizeof (double));
        select = ( a_int* ) malloc (ncv * sizeof (a_int));  // since HOWMNY = 'A', only used as workspace here

        a_int iparam[11], ipntr[11];
        iparam[0] = 1;                                      // ishift
        // iparam[2] = 10 * N;                                 // on input: maxit; on output: actual iteration
        iparam[2] = maxit;
        iparam[3] = 1;  // NB, only 1 allowed
        iparam[6] = 3;  // mode

        // char bmat[] = "G";                                  // generalized eigenvalue problem A*x = lambda*B*x
        // const char *bmat = "G";  // generalized eigenvalue problem A*x = lambda*B*x
        const char* bmat = convert_to_char (imp->opt.bmat);  // generalized eigenvalue problem A*x = lambda*B*x

        /*
c\Remarks
c  1. The converged Ritz values are always returned in ascending
c     algebraic order.  The computed Ritz values are approximate
c     eigenvalues of OP.  The selection of WHICH should be made
c     with this in mind when Mode = 3,4,5.  After convergence,
c     approximate eigenvalues of the original problem may be obtained
c     with the ARPACK subroutine dseupd .
        */
        // char which[]  = "LM";                               // compute the NEV largest (in magnitude) eigenvalues
        // const char* which    = imp->opt.which.c_str ();
        // char        howmny[] = "A";
        const char* which  = imp->opt.which.c_str ();
        const char* howmny = "A";

        // consider sigma = 0 for now.
#if defined(SG_USE_MKL)
        auto ls = getLinearSolver<Complex_t> ("CPardiso<Complex_t>");
#else
        auto ls = getLinearSolver<Complex_t> ("TPetscKsp<Complex_t>");
#endif

        if (sigma == 0.0)  // consider sigma = 0
        {
            ls->compute (K);
        }
        else  // consider sigma != 0
        {
            Complex_t      _sigma = sigma;
            ComplexSpMat_t Ksigma = K - _sigma * M;
            ls->compute (Ksigma);
        }

        a_int info = 0, ido = 0;

        // {
        //     // If INFO .EQ. 0, a randomly initial residual vector is used.
        //     // If INFO .NE. 0, RESID contains the initial residual vector, possibly from a previous run.
        //     info     = 1;
        //     auto eps = std::numeric_limits<double>::epsilon ();
        //     for (a_int i = 0; i < N; ++ i)
        //         resid[i] = eps;
        // }

        do
        {
            pznaupd_c (comm, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork, &info);
            // znaupd_c (&ido, convert_to_char(bmat), N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork,
            // &info);

            // solve Y = OP * X
            // dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
            if (ido == -1)
            {
                /*
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.
                */
                ComplexVec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x.set (i, ( Complex_t ) workd[ipntr[0] - 1 + i - row_start]);
                }
                x.assemble ();

                ComplexVec_t mx = M * x;

                ComplexVec_t y;
                ls->solve (mx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v = y[i];
                    memcpy (workd + ipntr[1] - 1 + i - row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido == 1)
            {
                /*
c          IDO =  1: compute  Y = OP * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    In mode 3,4 and 5, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.
                */
                ComplexVec_t mx (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    mx[i] = ( Complex_t ) workd[ipntr[2] - 1 + i - row_start];
                }
                mx.assemble ();

                ComplexVec_t y;
                ls->solve (mx, y);

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v = y[i];
                    memcpy (workd + ipntr[1] - 1 + i - row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido == 2)
            {
                /*
c          IDO =  2: compute  Y = B * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
                */
                ComplexVec_t x (N);
                for (int i = row_start; i < row_end; ++i)
                {
                    x[i] = ( Complex_t ) workd[ipntr[0] - 1 + i - row_start];
                }
                x.assemble ();

                ComplexVec_t mx = M * x;

                for (int i = row_start; i < row_end; ++i)
                {
                    Complex_t v = mx[i];
                    memcpy (workd + ipntr[1] - 1 + i - row_start, &v, sizeof (a_dcomplex));
                }
            }
            else if (ido != 99)
            {
                printf ("Error in znaupd: unexpected ido %d\n", ido);
                ret = 1;
                goto EXIT;
            }

        } while (ido != 99);

        printf ("znaupd: number of Arnoldi update iterations %d/%d\n", iparam[2], maxit);
        printf ("znaupd: NUMOP %d, NUMOPB %d, NUMREO %d\n", iparam[8], iparam[9], iparam[10]);

        // check info and number of ev found by arpack.
        if (info < 0 || iparam[4] < nev)
        {
            printf ("Error in znaupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
            ret = 1;
            goto EXIT;
        }

        pzneupd_c (comm,
                   rvec,
                   howmny,
                   select,
                   d,
                   z,
                   ldz,
                   sigma,
                   workev,
                   bmat,
                   N_local,
                   which,
                   nev,
                   tol,
                   resid,
                   ncv,
                   V,
                   ldv,
                   iparam,
                   ipntr,
                   workd,
                   workl,
                   lworkl,
                   rwork,
                   &info);

        if (info < 0)
        {
            printf ("Error in znaupd: info %d\n", info);
            ret = 1;
            goto EXIT;
        }

        // store the eigenvectors and eigenvalues
        eigenvalues.resize (nev);
        eigenvectors.resize (nev);

        for (int j = 0; j < nev; ++j)
        {
            eigenvectors[j].resize (N);

            for (int i = row_start; i < row_end; ++i)
            {
                eigenvectors[j].set (i, ( Complex_t ) z[j * N_local + i - row_start]);
            }

            eigenvectors[j].assemble ();

            eigenvalues[j] = ( Complex_t ) d[j];
        }

    EXIT:
        FREE_MEMORY (resid);
        FREE_MEMORY (V);
        FREE_MEMORY (z);
        FREE_MEMORY (d);
        FREE_MEMORY (workd);
        FREE_MEMORY (workl);
        FREE_MEMORY (select);

        return ret;
    }
}  // namespace

template <>
void ALGEBRA_EXPORT PArpack<Complex_t>::solve (const SparseMatrixType&  A,
                                               const SparseMatrixType&  B,
                                               std::vector<Complex_t>&  eigenvalues,
                                               std::vector<VectorType>& eigenvectors)
{
    imp->info = 0;
    imp->n    = A.rows ();
    // imp->nev   = 1;
    // imp->opt.sigma = 0.0;

    // checking the options...
    // see octave/liboctave/numeric/eigs-base.cc
    if (imp->opt.nev < 1 || imp->opt.nev > imp->n - 2)
    {
        printf ("Invalid number of eigenvalues to extract (must be 0 < k < n - 1) \n");
        return;
    }

    // [eigenvectors, eigenvalues]
    // eigenvectors: E(:. 0:nev)
    // eigenvalues: E(0:nev, nev)
    // ComplexMat_t E (imp->n, imp->opt.nev + 1);
    if (imp->opt.mode == arpack::Mode::MODE_2)
    {
        imp->info = eig2_z (imp, A, B, eigenvalues, eigenvectors);
    }
    else if (imp->opt.mode == arpack::Mode::MODE_3)
    {
        imp->info = eig3_z (imp, A, B, eigenvalues, eigenvectors);
    }

    return;
}